#' 所有结果一键分析
#'
#' @param dat harmonize后的dat
#' @param method_list 默认c("mr_wald_ratio","mr_egger_regression","mr_weighted_median"
#' ,"mr_ivw", "mr_simple_mode", "mr_weighted_mode")，可以用TwoSampleMR::mr_method_list()查看
#' 更多的方法
#' @param workers 多项目同时运行，加速分析，默认为1。
#' @param run_mr_presso 跑不跑mr_presso，默认为TRUE
#' @param NbDistribution 默认3000，可以1000,10000
#'
#' @return res
#' @export
#'
#' @examples
#'
#' \dontrun{
#' res <- U5_mr(dat)
#' }
#'

U5_mr <- function(dat,
                  method_list = c("mr_wald_ratio","mr_egger_regression","mr_weighted_median"
                                  ,"mr_ivw", "mr_simple_mode", "mr_weighted_mode"),
                  workers=1,
                  run_mr_presso=TRUE,
                  NbDistribution=3000)
{
  dat <- subset(dat, mr_keep) %>%
    split(list( .$id.exposure , .$id.outcome) )

  dat <- Filter(function(x){ nrow(x) > 0} , dat)

  future::plan( future::multisession, workers = workers )

  message( "\n\n分析基本方法中......" )
  mr_tab <- dat  %>%
    furrr::future_map_dfr( mr_to_wide,method_list=method_list ,.options = furrr::furrr_options(seed = TRUE),.progress = TRUE
    )

  if(run_mr_presso){
  message( "\n\n分析MR_PRESSO中......" )

  presso <- dat %>%
    furrr::future_map_dfr( add_presso,NbDistribution=NbDistribution ,.options = furrr::furrr_options(seed = TRUE),.progress = TRUE
    )
  mr_tab <- dplyr::full_join(mr_tab,presso, by = c('id.exposure','id.outcome'))
  }


  message( "计算F值，R2中......" )
  message("使用2*eaf*(1-eaf)*beta^2/(2*eaf*(1-eaf)*beta^2+ 2*eaf*(1-eaf)*effective_n*se^2  )的公式计算R2")
  message("使用beta^2/se^2公式计算单个F(singleF)值")


  Fst <- dat %>%
    furrr::future_map_dfr( add_F ,.options = furrr::furrr_options(seed = TRUE),.progress = TRUE
    )
  mr_tab <- dplyr::full_join(mr_tab,Fst, by = c('id.exposure','id.outcome'))


  message( "\n\n计算异质性，多效性中......" )

  het_pleio <- dat %>%
    furrr::future_map_dfr( add_het_pleio ,.options = furrr::furrr_options(seed = TRUE),.progress = TRUE
    )
  mr_tab <- dplyr::full_join(mr_tab,het_pleio, by = c('id.exposure','id.outcome'))




  mr_tab <- mr_tab %>%
    dplyr::mutate_at(dplyr::vars(dplyr::starts_with( c("b_","se_","nsnp_","pval_","lo_ci_","up_ci_",
                                                       "or_","or_lci_95_","or_uci_95_",
                                                       "total_R2","Fstat","mF","minF",
                                                       "maxF","IVW_Q","Egger_Q","egger_intercept"))), as.numeric)

  message( "\n\n计算统计效能中......" )

  samplesize <- dat %>%
    furrr::future_map_dfr( add_samplesize ,.options = furrr::furrr_options(seed = TRUE),.progress = TRUE
    )
  mr_tab <- dplyr::full_join(mr_tab, samplesize, by = c('id.exposure','id.outcome'))


  mr_tab$b1 <- abs( ifelse( !is.na(mr_tab$`b_Inverse variance weighted`),
                            mr_tab$`b_Inverse variance weighted`, mr_tab$`b_Wald ratio` ) )



  for (i in 1:nrow(mr_tab) ) {

  if(mr_tab$type.outcome[i]=="Binary" ){

    mr_tab$power[i] <-pnorm(sqrt(mr_tab$samplesize.outcome[i]*mr_tab$total_R2[i]*(mr_tab$ratio[i]/(1+mr_tab$ratio[i]))*(1/(1+mr_tab$ratio[i])))*mr_tab$b1[i]-qnorm(1-0.05/2))

    mr_tab$power_b[i]<-((qnorm(0.8) + qnorm(1-0.05/2))/sqrt(mr_tab$samplesize.outcome[i]*mr_tab$total_R2[i]*(mr_tab$ratio[i]/(1+mr_tab$ratio[i]))*(1/(1+mr_tab$ratio[i]))))
  }else{
    mr_tab$power[i] <-pnorm(sqrt(mr_tab$samplesize.outcome[i]*mr_tab$total_R2[i])*mr_tab$b1[i]-qnorm(1-0.05/2))

    mr_tab$power_b[i]<-((qnorm(0.8) + qnorm(1-0.05/2))/sqrt( mr_tab$samplesize.outcome[i] * mr_tab$total_R2[i] ) )
  }
  }


  message( "计算统计效能power文献为https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4052137/" )

  message( "计算统计效能power上下限文献为https://www.nature.com/articles/s41593-022-01174-7#Sec13" )

  mr_tab$power_lower_beta <- - mr_tab$power_b
  mr_tab$power_upper_beta <- mr_tab$power_b
  mr_tab$power_lower_or<- exp(- mr_tab$power_b)
  mr_tab$power_upper_or <- exp( mr_tab$power_b)

  mr_tab$required_beta_to_meet_0.8power <- sprintf(" ≤ %.3f or ≥ %.3f ", mr_tab$power_lower_beta, mr_tab$power_upper_beta)
  mr_tab$required_OR_to_meet_0.8power <- sprintf(" ≤ %.3f or ≥ %.3f ", mr_tab$power_lower_or, mr_tab$power_upper_or)


  return(mr_tab)
}


mr_to_wide <- function(dat,method_list)
{
  suppressPackageStartupMessages(require(TwoSampleMR,quietly = TRUE))
  x <- subset(dat, mr_keep)
  x <- x %>% dplyr::distinct(SNP,.keep_all = TRUE)
  if(nrow(x) == 0)
  {
    message("No SNPs available for MR analysis of '", dat$id.exposure[1], "' on '", dat$id.outcome[1], "'")
    return(NULL)
  } else {
    # message("\n Analysing '", dat$id.exposure[1], "' on '", dat$id.outcome[1], "'")
  }

  res <- lapply(method_list, function(meth)
  {
    get(meth)(x$beta.exposure, x$beta.outcome, x$se.exposure, x$se.outcome, TwoSampleMR::default_parameters())
  })
  methl <- TwoSampleMR::mr_method_list()
  mr_tab <- data.frame(
    outcome = x$outcome[1],
    exposure = x$exposure[1],
    method = methl$name[match(method_list, methl$obj)],
    nsnp = sapply(res, function(x) x$nsnp),
    b = sapply(res, function(x) x$b),
    se = sapply(res, function(x) x$se),
    pval = sapply(res, function(x) x$pval)
  )

  mr_tab <-TwoSampleMR::generate_odds_ratios(mr_tab)

  mr_tab$`Beta (95% CI)` <- ifelse(is.na(mr_tab$b), NA,
                                   sprintf("%.3f (%.3f to %.3f)",
                                           mr_tab$b, mr_tab$lo_ci, mr_tab$up_ci))


  mr_tab$`OR (95% CI)` <- ifelse(is.na(mr_tab$or), NA,
                                 sprintf("%.3f (%.3f to %.3f)",
                                         mr_tab$or, mr_tab$or_lci95, mr_tab$or_uci95))
  mr_tab1<-c(t(mr_tab))
  mr_tab2<-as.data.frame(mr_tab1)
  mr_tab3<-as.data.frame(t(mr_tab2))

  name<-c("outcome","exposure","method","nsnp","b","se","pval","lo_ci","up_ci","or","or_lci95","or_uci95","Beta (95% CI)","OR (95% CI)")
  name_all<-c()
  for (i in 1:length(methl$name[match(method_list, methl$obj)])) { name2<- paste0(name,"_",(methl$name[match(method_list, methl$obj)])[i])
  name_all<-c(name_all,name2)
  }

  colnames(mr_tab3)<-name_all

  mr_tab3$exposure<- dat$exposure[1]
  mr_tab3$outcome<- dat$outcome[1]
  mr_tab3$id.exposure<- dat$id.exposure[1]
  mr_tab3$id.outcome<- dat$id.outcome[1]
  mr_tab3$nSNP<- length(dat$exposure)

  return(mr_tab3)
}


add_presso<- function(dat,NbDistribution){
  dat <- subset(dat, mr_keep)
  dat <- dat %>% dplyr::distinct(SNP,.keep_all = TRUE)
  if(length(dat$SNP) >3){

    PRESSO <-MRPRESSO::mr_presso(BetaOutcome ="beta.outcome", BetaExposure = "beta.exposure", SdOutcome ="se.outcome", SdExposure = "se.exposure",
                                 OUTLIERtest = TRUE,DISTORTIONtest = TRUE, data = dat, NbDistribution = NbDistribution,
                                 SignifThreshold = 0.05)
    exposure<-dat$exposure[1]
    outcome<-dat$outcome[1]
    method<-"MR_PRESSO"
    if( is.na(PRESSO[["Main MR results"]][["Sd"]][2]) ){n_outliers<-0
    }else{n_outliers<-length(dat$SNP[PRESSO[["MR-PRESSO results"]][["Distortion Test"]][["Outliers Indices"]]])}
    nsnp<-length(dat$SNP)-n_outliers
    b<-PRESSO[["Main MR results"]][["Causal Estimate"]][2]
    se<-PRESSO[["Main MR results"]][["Sd"]][2]
    pval<-PRESSO[["Main MR results"]][["P-value"]][2]

    PRESSO_res<- as.data.frame(cbind(outcome,exposure,method,nsnp,b,se,pval))
    PRESSO_res$b<- as.numeric(PRESSO_res$b)
    PRESSO_res$se<- as.numeric(PRESSO_res$se)
    PRESSO_res$pval<- as.numeric(PRESSO_res$pval)
    PRESSO_or<- generate_odds_ratios(PRESSO_res)

    Global_Test_RSSobs <-  PRESSO[["MR-PRESSO results"]][["Global Test"]][["RSSobs"]]

    Global_Test_Pvalue <- as.character(PRESSO[["MR-PRESSO results"]][["Global Test"]][["Pvalue"]])

    if( !is.na(PRESSO_or$b[1]) ){

      PRESSO_or$`Beta (95% CI)` <- ifelse(is.na(PRESSO_or$b), NA,
                                          sprintf("%.3f (%.3f to %.3f)",
                                                  PRESSO_or$b, PRESSO_or$lo_ci, PRESSO_or$up_ci))


      PRESSO_or$`OR (95% CI)` <- ifelse(is.na(PRESSO_or$or), NA,
                                        sprintf("%.3f (%.3f to %.3f)",
                                                PRESSO_or$or, PRESSO_or$or_lci95, PRESSO_or$or_uci95))

      colnames(PRESSO_or)<-paste0(colnames(PRESSO_or),"_MR_PRESSO")

      outliers<-dat$SNP[PRESSO[["MR-PRESSO results"]][["Distortion Test"]][["Outliers Indices"]]]

      #把离群值放到一个格子里
      outliers<-paste(outliers,sep = ",",collapse=" ")
      PRESSO_Distortion_beta<-PRESSO[["MR-PRESSO results"]][["Distortion Test"]][["Distortion Coefficient"]][["beta.exposure"]]
      PRESSO_Distortion_Pvalue<- as.character(PRESSO[["MR-PRESSO results"]][["Distortion Test"]][["Pvalue"]])

      presso_OR<-cbind(PRESSO_or,Global_Test_RSSobs,Global_Test_Pvalue,outliers,PRESSO_Distortion_beta,PRESSO_Distortion_Pvalue)

    }else{
      colnames(PRESSO_or)<-paste0(colnames(PRESSO_or),"_MR_PRESSO")
      presso_OR<-cbind(PRESSO_or,Global_Test_RSSobs,Global_Test_Pvalue)
    }


  }else{

    message(dat$id.exposure[1],"到",dat$id.outcome[1],"的SNP个数少于4，无法计算PRESSO")

    exposure_MR_PRESSO<-dat$exposure[1]
    outcome_MR_PRESSO<-dat$outcome[1]
    method_MR_PRESSO<-"MR_PRESSO"
    nsnp_MR_PRESSO<-length(dat$SNP)

    presso_OR<-as.data.frame(cbind(outcome_MR_PRESSO,exposure_MR_PRESSO,method_MR_PRESSO,nsnp_MR_PRESSO))

  }
  presso_OR$id.exposure<-dat$id.exposure[1]
  presso_OR$id.outcome<-dat$id.outcome[1]
  return(presso_OR)

}


add_F <- function(dat){

  dat <- subset(dat, mr_keep)
  dat <- dat %>% dplyr::distinct(SNP,.keep_all = TRUE)
  if( ( "ncase.exposure" %in% colnames(dat)) & ( "ncontrol.exposure" %in% colnames(dat) ) ){
    dat$ncase.exposure<-as.numeric(dat$ncase.exposure)
    dat$ncontrol.exposure<-as.numeric(dat$ncontrol.exposure)
    dat$effective_n.exposure<- ifelse( !is.na(dat$ncase.exposure) & !is.na(dat$ncontrol.exposure),
                                       TwoSampleMR::effective_n(dat$ncase.exposure,dat$ncontrol.exposure ),
                                       dat$samplesize.exposure
    )
  }else{
    dat$effective_n.exposure<- dat$samplesize.exposure
  }

  if("eaf.exposure" %in% colnames(dat) ){


    dat$R2<-with(dat,
                 (2*eaf.exposure*(1-eaf.exposure)*beta.exposure*beta.exposure)/(2*eaf.exposure*(1-eaf.exposure)*beta.exposure*beta.exposure+2*eaf.exposure*(1-eaf.exposure)*effective_n.exposure*se.exposure*se.exposure)
    )
  }else if("maf.exposure" %in% colnames(dat)){
    dat$R2<-with(dat,
                 (2*maf.exposure*(1-maf.exposure)*beta.exposure*beta.exposure)/(2*maf.exposure*(1-maf.exposure)*beta.exposure*beta.exposure+2*maf.exposure*(1-maf.exposure)*effective_n.exposure*se.exposure*se.exposure)
    )
  }else{
    message("需要maf或eaf计算R2")
    dat$R2 <- NA
  }



  dat <- dat %>%
    dplyr::mutate(singleF = beta.exposure^2/(se.exposure^2) )

  mF = mean(dat$singleF)

  minF = min(dat$singleF)
  maxF = max(dat$singleF)

  #计算R2
  total_R2<-sum(dat$R2)

  #计算F值
  Fstat <- with(dat,
                (effective_n.exposure-length(SNP)-1)/(length(SNP)) * total_R2/(1-total_R2)
  )
  Fstat <-max(Fstat)


  F_R2<-as.data.frame(cbind(total_R2,Fstat,mF,minF,maxF))

  F_R2$id.exposure<-dat$id.exposure[1]
  F_R2$id.outcome<-dat$id.outcome[1]

  return(F_R2)
}


add_het_pleio <- function(dat){

  suppressPackageStartupMessages(require(TwoSampleMR,quietly = TRUE))

  dat <- subset(dat, mr_keep==TRUE)
  dat <- dat %>% dplyr::distinct(SNP,.keep_all = TRUE)

  if(length(dat$SNP)>2){

    het <- mr_heterogeneity(dat)

    het_IVW<-subset(het,method== "Inverse variance weighted")


    if(!identical(het,het_IVW )){

      het_IVW <- subset(het, method == "Inverse variance weighted")
      het_Egger <- subset(het, method == "MR Egger")
      het_out <- as.data.frame(cbind(het_IVW[, 5:8], het_Egger[,
                                                               5:8]))
      colnames(het_out) <- c("method_IVW", "IVW_Q", "IVW_Q_df",
                             "IVW_Q_pval", "method_Egger", "Egger_Q", "Egger_Q_df",
                             "Egger_Q_pval")

    }else{
      het_IVW <- het
      method_Egger <- "MR Egger"
      het_out <- as.data.frame(cbind(het_IVW[, 5:8], method_Egger))
      colnames(het_out) <- c("method_IVW", "IVW_Q", "IVW_Q_df",
                             "IVW_Q_pval", "method_Egger")
    }


  }else if(length(dat$SNP)>1){
    het <- mr_heterogeneity(dat)

    het_IVW<-het
    method_Egger<-"MR Egger"


    het_out<-as.data.frame(cbind(het_IVW[,5:8],method_Egger))
    colnames(het_out)<-c("method_IVW","IVW_Q","IVW_Q_df","IVW_Q_pval","method_Egger")

  } else {

    method_IVW<-"Inverse variance weighted"

    method_Egger<-"MR Egger"


    het_out<-as.data.frame(cbind(method_IVW,method_Egger))
  }


  #多效性
  if(length(dat$SNP)>2){

    pleio <- mr_pleiotropy_test(dat)

    pleio<-pleio[,5:7]
    egger_pleio<-as.data.frame(pleio)
    colnames(egger_pleio)<-c("egger_intercept","egger_intercept_se","egger_intercept_pval")

  }else{
    egger_pleio<-c()

  }

  out2<-as.data.frame(c(het_out,egger_pleio))


  out2$id.exposure<-dat$id.exposure[1]
  out2$id.outcome<-dat$id.outcome[1]

  return(out2)
}


add_samplesize<-function (dat)
{
  dat <- subset(dat, mr_keep)
  dat <- dat %>% dplyr::distinct(SNP, .keep_all = TRUE)
  out <- data.frame(id.exposure = dat$id.exposure[1],
                    id.outcome = dat$id.outcome[1])

  out$samplesize.exposure <- dat$samplesize.exposure[1]
  if( "ncase.exposure" %in% colnames(dat)){ out$ncase.exposure <- dat$ncase.exposure[1] }
  if( "ncontrol.exposure" %in% colnames(dat)){ out$ncontrol.exposure <- dat$ncontrol.exposure[1]}
  out$type.exposure <- dat$type.exposure[1]

  out$samplesize.outcome <- dat$samplesize.outcome[1]
  if( "ncase.outcome" %in% colnames(dat)){ out$ncase.outcome <- dat$ncase.outcome[1] }
  if( "ncontrol.outcome" %in% colnames(dat)){ out$ncontrol.outcome <- dat$ncontrol.outcome[1]}
  out$type.outcome <- dat$type.outcome[1]
  out$ratio <-dat$ratio[1]
  return(out)
}





